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Abstract 

An approximate method based on adiabatic time dependent density functional the- 

ory (TDDFT) is presented, that ahows for the description of the electron dynamics in 

nanoscale junctions under arbitrary time dependent external potentials. In this scheme, 

H the density matrix of the device region is propagated according to the Liouville-von Neu- 

^ mann equation. The semi-infinite leads give rise to dissipative terms in the equation of 

motion which are calculated from first principles in the wide band limit. In contrast to 
Jr^ earlier ab-initio implementations of this formalism, the Hamiltonian is here approximated 

by a density expansion in the spirit of the density functional based tight-binding (DFTB) 
^ method without introducing empirical parameters. Results are presented for two proto- 

^ typical molecular devices and compared to calculations at the full TDDFT level. The 

issue of non-existence of a steady state under certain conditions is also briefly touched 
S on. 

o 

O !• Introduction 

^ Molecular electronics - the use of single molecules as functional entities in electronic 

^ devices - is often heralded as a replacement for the conventional CMOS technology which 

^\ faces physical limits in further miniaturization. Even if the industrial application of this 

CN| technology seems to be out of reach in the moment, research on the electronic transport 

through individual molecules has already revealed a variety of interesting quantum phe- 
nomena. To these belong strong electron correlation giving rise to Coulomb blockade and 
Kondo physics, elect ron-phonon interactions responsible for device heating and measur- 
able in inelastic tunneling spectra, as well as electron-photon interactions for molecular 
I photo-switches [l]. While a qualitative understanding of these complex processes has 

• • been obtained in many cases, quantitative agreement between first-principles theory and 

, ^ experiment is still unsatisfactory for the seemingly simple case of ballistic steady state 

transport [2^. 
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The majority of such studies for reahstic devices are based on the nonequihbrium 
Green's function (NEGF) technique [3 , which reduces to the traditional Landauer trans- 
mission formahsm [HE] for coherent transport. The device Green's function is then con- 
structed from the Kohn-Sham (KS) single-particle Hamiltonian of ground state Density 
Functional Theory (DFT) including the effects of the leads through suitable self-energies. 
For finite bias, more sophisticated approaches determine the actual potential profile on 
the molecule by solving the Poisson equation in order to achieve a self-consistent treat- 
ment of Hamiltonian and Green's function [1 . 

In recent years, also time-dependent Density Functional Theory (TDDFT) has be- 
come popular in the context of molecular electronics j6HT3]. The reason for this is 
twofold. First, TDDFT goes beyond the ad-hoc application of the DFT Hamiltonian 
for non-equilibrium systems and provides a more rigorous theoretical foundation [M]. 
In fact, calculations at this level provide corrections to the conventional DFT-NEGF 
results for steady state currents, provided non-local exchange-correlation functionals are 
employed [T5HIH]. Second, TDDFT allows for the description of dynamical processes like 
the switching behavior of molecular devices and alternating currents. The formalism may 
also be easily extended to cover the interaction with light, which paves the way to study 
photo- assisted and photo-suppressed transport, the fiuorescence of contacted molecules 
or atomistic opto-electronic devices like photoswitches. 

So far, several implementations of TDDFT for open systems have been suggested, 
which may be classified according to the way the transport boundary conditions are 
treated. The simplest approach is to approximate the mesoscopically large leads by 
atomic clusters of finite size O |9l [TTJ [13]. The Kohn-Sham orbitals are then time- 
propagated under the infiuence of a potential difference between left and right leads and 
changes of molecular charges in the contacts allow one to quantify the current. This 
approach has the advantage that TDDFT algorithms for isolated systems can be applied 
to the transport scenario without large changes. On the downside, a steady-state current 
is only transiently reached due to charge accumulation. Also, unphysical oscillations of 
the current never fully die out, because of the discrete level spectrum of the contacts. 

Another path was pursued by Burke and co-workers [H [19] . Here a ring-topology of 
the electronic circuit is used to avoid the conceptual difficulty of having two different 
chemical potentials as in the Landauer picture. A constant electric field on the ring is 
used to generate a current and coupling to a heat bath within a master equation approach 
allows for the establishment of a steady-state. 

A third approach is given by a separation of the conductor into a finite device region 
and two semi-infinite leads, very similar to the treatment within the standard DFT- 
NEGF approach (see Fig. [l^). The time-dependent Kohn-Sham equations need to be 
solved for the central region only and the effect of the leads is exactly accounted for by 
properly defined self-energies. Several model studies along these lines are documented 
in the literature [TOl [201422] . Finally, a closely related path is followed by Chen and co- 
workers [12 , 23]l24l. Instead of propagating Kohn-Sham states, the equations of motion 
are written in terms of the reduced density matrix for the device region. This fact allows 
one to obtain the initial conditions from a straightforward application of equilibrium 
Green's function theory. 

For all of the methods mentioned above, simulations on realistic molecular devices 
which take the atomistic structure of both molecule and leads into account are numer- 
ically demanding. This is because the electron motion has to be fully resolved, leading 
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Figure 1: (a) Schematic device setup with D denoting the central device region. The left (L) and right 
(R) leads are assumed to be semi-infinite and in thermal equilibrium, (b) Seven-membered carbon chain 
between Al nanowires in (001) direction. Frontier carbon atoms in hollow position of the Al surface at 
1.0 A distance. Bond lengths (A): C-C = 1.32, Al-Al = 2.86. (c) 1,4-benezenediol molecule inbetween 
Al nanowires. 



to time steps in the sub-fs range. Although only the central part of the device is out of 
equilibrium, this region has to include several layers of the lead material in addition to 
the functional molecule in question, which raises the dimension of the Hamiltonian. Only 
in this way a smooth transition of potentials at the device-lead interface is guaranteed. 
Moreover, comprehensive studies require numerous simulations involving the variation 
of different parameters, like different molecular binding configurations, temporal bias 
profiles or light frequency and amplitude in case of photo-induced processes. Hence, 
a computationally efficient, but still predictive method based on TDDFT seems to be 
desirable. 

In this contribution we present efforts in this direction and discuss an implementation 
of the density functional based tight-binding method (DFTB) [25H27] for open systems. 
The ground state DFTB approach is characterized by several additional approximations 
beyond a particular exchange-correlation functional and has already been generalized 
to the time domain (TD-DFTB, |28ti3Q]) in the spirit of TDDFT. Here we follow the 
mentioned density matrix approach of Chen and co-workers to arrive at a fast simulation 
method for open systems driven by time dependent bias potentials. In Sec. [2] the basic 
equations of motion are reviewed, while Sec. [3] discusses the approximations underlying 
the DFTB scheme and the necessary adaptions in the present context. Several case 
studies are presented in Sec. [4] followed by a brief summary and outlook in Sec. [5] 

2. Equation of motion for the reduced density matrix 

2. 1 . General formulas 

We consider a system as depicted in Fig. [l^ and aim at evaluating the current through 
the central device region (D) which is driven by a time-dependent bias potential V{t). 
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The left (L) and right (R) lead are thought to extend from the device region to infin- 
ity and feature a continuous density of states to ensure proper dissipation. In second 
quantization, the electronic Hamiltonian of the whole system can be written as 

mn a,k a,k;m 

Here, the first term describes the device region, where dl^{dm) is the creation (annihi- 
lation) operator for an electron in the atomic orbital m, and Hmn stands for the corre- 
sponding Hamiltonian in this basis. The second term describes the contacts (a G {L, 
R}), where cj^ f.{ca,k) is the creation (annihilation) operator for an electron in the single- 
particle state k in the ath contact, while ea,k is the corresponding level energy. The 
third term describes the coupling between device region and contacts and features the 
hopping matrix T (matrices will be denoted by bold face letters). Although the Hamilto- 
nian in Eq. ([T]) is formally valid only for non-interacting particles, the yet to be specified 
parameters may effectively account for electron-electron interactions, for example in a 
DFT context. Moreover, a self consistent determination of the device Hamiltonian and 
hopping matrix will result in a time dependence of these quantities for varying external 
bias. 

Next, we define the lesser Green's function [31 of the device region as 

G<^{t,t') = i{dl{t')d^{t)), (2) 
and apply the equation of motion for Heisenberg operators to arrive at 

i^G< {t, t') = [H, G< {t, t')] + E it, t'), (3) 

a 

where Q^{t^t') is given in terms of the coupling matrix and the lead-device lesser 
Green's function 

Glu;nM=i{di{t')c^,k{t)) (4) 

as follows 

k 

Following the standard NEGF technique [ST, the term Q^it) can be further rewritten 
in a form that incorporates only matrices with the dimension of the device region: 

/ + 00 
dr[S^(t,T)G<(r,i') + S<(t,r)G"(r,0] 
-CO 



CO 
+ CO 



dr[G< {t, r)K{r, t') + G^{t, t)S< (r, t')]. (6) 



Here, G^*^^^ are the retarded (advanced) Green's functions of the device region, which 
satisfy the Dyson equation 

(i|-ifW)G'-(°)(t,t') 

/ + 00 
dri:i^'^Ht,r)G^^'^\r,t')=5{t-t'), (7) 
-CO 



while 5]^'*^'^ are the retarded (advanced, lesser) self-energies due to the coupling to the 
ath contact. They are given as 

Sr><(t,r) = Tt (t)g^.».<(t,r)r«(r), (8) 

with g^^''^ denoting the Green's functions of the ath contact. 

Equations ( 3|6|7|8 ) form a closed set of equations to describe the dynamical properties 



of non-equilibrium systems and will be the theoretical basis of our method. Instead of 
working with the full lesser Green's function, a significant simplification is achieved by 
turning to the reduced density matrix for the device region {cr{t)) as the key quantity, 
which is dependent on a single time argument only. The quantity a{t) is directly related 
to the lesser Green's function G^{t^t^) 

a{t) = (9) 

Thus, the equation of motion for a{t) can be obtained from equation ([3| as 

i^^a{t) = [H{t),a{t)]-i^Q^{t). (10) 



For vanishing coupling (i.e. T = 0), Q vanishes and equation (10) reduces to the con- 
ventional quantum Liouville equation for closed systems. The quantity describes 
dissipation effects due to the presence of the semi-infinite contacts and gives rise to finite 
lifetimes of states located in the central region. It is also noted that the trace of the 
matrix gives the particle current through the interface Sa (see Fig. [l^) in contact a 

m 

Ut) = -Tr[QJt)]. (11) 

2.2. Wide-hand approximation for the dissipation term 

Direct application of Eq. ^ to evaluate the dissipation term turns out to be 
computationally demanding. This is due to the necessity to store and integrate over 
two-time quantities like the self-energies. A common solution is to invoke the wide band 
approximation (WBA), where both the lead density of states and the device- lead coupling 
are assumed to be smooth and not strongly dependent on energy. This is a reasonable 
approximation for simple metals in the linear response regime but fails for larger applied 
bias and/or leads with a more complex density of states, like for example nanotubes. 
The WBA is however not without alternative in the present context. Recently, Zheng 
and co-workers presented a hierarchical equation of motion approach, which goes beyond 
the WBA at a tolerable increase in computational effort [24J. 

Within the WBA, the self-energies X)^^^ become local in time and can be written as 

m 

l^TM = (r„TiA„)<5(i-0- (12) 

Here, the matrices Fq, and Aq, denote the hermitian and anti-hermitian part of 
5]^(£^), respectivel}j^ and are evaluated at the common Fermi energy {Ep) of the un- 
biased system. Assuming that the non-interacting electrons in contact a are in local 



■•^Note that the definition of Fq, and Aq, is often interchanged in the hterature. We keep with the 
notation of reference |12| to ease comparison. 
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equilibrium characterized by a Fermi-Dirac distribution function /"(e) with chemical 
potential /io and assuming further that the bias potential Va{t) leads to a rigid shift of 
the lead energy levels, one obtains 

S<(t,0 = ^A«exp|i^*dTT/«(T)|y^°°der(e)e--(*-*'). (13) 

With the expressions ( 12|13 ) above, the dissipation term Qa{t) can be simplified to 

[E] 

Q^{t) = i[T^,a{t)] + {A„ a{t)} + K^{t). (14) 

Here, the first term represents the renormalization of the device energy levels due to the 
ath contact and involves a commutator; the second term involving an anti-commutator 
describes the level broadening, while the hermitian matrix Ka{t) 

2i^^ _ r° dee' 



/ N 2i ^ , , dee''"^ 



-oo 

2i ... ,,,, de 



. J- J - ^"^^^^'"J (e - V.it)) I - Hit) - E- ^- + (1^) 



incorporates the history of the applied bias via the propagator Ua{t) 
U^{t) = exp j\H{T) + S'- - IV^{T)]dT^ 



(16) 



The term 5]^ is a summation of the retarded self-energy functions for both contacts. 

Equations ([Toj [l4j [Tsj [T6| provide an approximate approach to calculate the transient 
non-equilibrium density matrix for open systems. Zheng et. al. have combined the scheme 
described above with TDDFT to simulate the time-dependent transport in nanoscale 
systems |12][23l. Next, we will describe how the density matrix cr(t), Hamiltonian H{t)^ 
and dissipation term Q^it) are constructed in the framework of the density- functional 
based tight-binding (DFTB) method, which aims to provide a more efficient but still 
accurate method for real material simulations. 



3. Density Functional based Tight-Binding Method 

3.1. DFTB for closed systems 

As an approximate DFT method, DFTB was initially developed to study the ground 
state of finite or periodic systems. Details on the derivation and performance for these 
class of materials may be found in several reviews [32-34 . As in conventional tight- 
binding models, the total energy of DFTB consists of a band structure part £^bs and a 
repulsive energy part £^rep7 

E^^^^=E^s+E,,^. (17) 

Since we will not consider structural changes during electron transport in the present 
work, only the band structure term is relevant for the further discussion. It takes the 
form 
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occ ^ 
^BS = ^{^^\Ho\^^) + - ^ 7c./3 Ag« A^;^ , (18) 



2 

a, (3 



which is obtained through an expansion of the DFT total energy around a reference 
density n^{r) up to second order [27 . The latter is taken to be a superposition of the 
neutral densities of all atoms comprising the system in question (n^(r) = Xla^al^))- 



The first term in Eq. ( 18) involves a sum over expectation values of the DFT Hamiltonian 
evaluated at n^: 

Hoir) = ~ + Ve.t{r) + J dr'^^^^ + V^c[no], (19) 

with the kinetic energy operator, the external potential, the Hartree potential generated 
by no and the DFT exchange-correlation potential. 



The last quantity in Eq. (18) represents the second-order term in the mentioned 
density expansion and accounts for charge transfer effects. Here, Aq^ is defined as the 
difference of the Mulliken charges between charged and neutral atom a, i.e. Aq^ = 
Qa — Qa- The matrix j^f^ is a measure of the electron-electron interaction and decays 
like l/\Ra — R(3\ for large distances between atoms at and Rf^. For the on-site case, 
the Hubbard-like parameter Ua = ^aa is taken from full atomic DFT calculations and 
represents the chemical hardness of the respective element. An interpolation formula 
between these limiting cases can be analytically constructed t>y invoking a monopole 
approximation for the density fluctuations Sria, which are assumed to take the form 

Sna{r) = ^^e-^"l^-^"l, with / dr Sna{r) = Aq^. (20) 

OTT J 



Here, the effective decay constant Tq, = is determined by the constraint Ua = 7( 



In the practical implementation of DFTB, the Kohn-Sham orbitals 1^^) are expanded 
in a set of localized atomic orbitals (AO) that are again obtained from a full atomic 

DFT calculation 

\^i) = Y.c,^\^^,). (21) 

Usually only occupied AO are used in the molecular expansion, but this is not a necessary 
requirement of the model. The atomic DFT calculations are converged with respect to the 
basis set size, so that each AO is constructed from a superposition of many Slater-type 
orbitals. Furthermore, an additional confinement potential in the atomic calculations 
leads to a suppression of the long-range tail of the AO, making them suitable expansion 
functions for the molecular case [26]. The minimal basis set of DFTB is expected to 
be more accurate than conventional contracted basis sets of the same size, e.g. Pople's 
ST0-3G. 



In a two-center approximation, matrix elements of the DFT Hamiltonian in Eq. (19) 
{H^j^) are numerically computed and stored in Slater-Koster [35 tables as function of 
the inter-element distance. The Mulliken charges used to estimate the second-order term 
in the density expansion are given in terms of the molecular orbital coefficients (q^) and 
the overlap of two atomic orbitals {Sjj^u = (0^|0j,)): 
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^ OCC 

~ 2 ^ ^ + ^lyCi^Sy^). (22) 

The trace is performed over all atomic orbitals on site a and cr denotes the density matrix 
in the AO basis. Thus the total energy of DFTB will be a function of the expansion 
coefficients {ci^^c*^} in the given basis set. The ground state energy can be obtained 
from the variational principle 

S{Ebs - ^iC*^Ci,S^,)/5c*^ = 0, (23) 

where the are Lagrange multipliers. 

The solution of equation ( [23| ) gives rise to the secular equation 

Y^^H^, - e,S^,)c,, = 0, (24) 
where the Hamiltonian matrix element i^^^^ for atom pair (a, /3) reads 

Hf^iy = + ^^M^ X^(7c.7 + 7^7)^^7 I^^Oi.v e p. (25) 



Equation (24) can be viewed as an approximate version of the Kohn-Sham equation 



in DFT. The second term in (25) accounts for changes of the electron density with 
respect to the reference (n^) and needs to be determined self-consistently via repeated 
evaluation of Eqs. ( [24| [22| and[25|), very similar to self consistent field approaches. The 
computational efficiency of the DFTB method grounds on two facts. First, the minimal 
but still accurate basis set is smaller then in conventional DFT approaches and second, 
the construction of the Hamiltonian requires little CPU time due the precomputation of 
all numerical integrals. 

3.2. DFTB for open systems 

For open systems, the concept of a total energy is not well-defined, and the varia- 
tional principle and secular equation do not hold. Thus, the DFTB method needs to be 
generalized in oder to deal with a semi-infinite system out of equilibrium. Recalling the 
equations ( 1Q|14|15|16[ ) in Sec. [5] which describe the transient non-equilibrium dynamics 



in open systems, the key ingredients which determine the reduced density matrix a(t) are 
the Hamiltonian H{t) and self-energy functions Xlc^. Here we present how to construct 
these matrices in the framework of DFTB. 

3.2.1. Device Hamiltonian 

In the spirit of DFTB, the Hamiltonian should take the general form 

H^, = HlM+5V^A{^qc,}] (26) 



where the first term is the zero-order DFTB Hamiltonian, while the second term includes 
the Hartree potential SVh and exchange-correlation potential SVxc due to the difference 
density 5n = Sua. 

For closed systems, where the electrostatic potential can be assumed to tend to zero 
at infinity, the Hartree potential SVnir) is given by 

6VH{r) = I dr'p^^. (27) 



Using Eq. (20), this results in the 7- matrix dependent second-oder term of the Hamilto- 
nian in Eq. (25 ). For an electronic device like the one depicted in Fig.[l^, the electrostatic 
potential in the left and right contact is set by the applied bias. The potential in the 
central region must therefore be obtained from the Poisson equation 

V'^dVHir) = -47r(5n(r) (28) 

under these boundary conditions. In our approach, the device density matrix (T{t) is 



used to compute time-dependent Mulliken charges according to Eq. (22), which are used 



to set up the difference density 8n{r) via Eq. (20). Taking Tq, = WUa as above, some 



exchange and correlation effects are taken into account at this point, since the Hubbard- 
like parameters Ua stem from a DFT calculation and are considerably smaller than 
unscreened Hartree-only parameters [36] . 

The Poisson equation is then discretized on a grid and solved in real space for a 
sufficiently large box enclosing the device region. Note that an inherent assumption of 
such an approach is that the potential at the device-lead interface equals the potential 
in the contact interior. In other words, the device region should include enough layers of 
the contact material to ensure efficient screening. 

Similar to the electron density, 5VH{r) is projected on the atomic sites through 



5Vc, = \ [ drSVH{r)e-^-\'^-^-\ (29) 

^ -/box 



with normalization factor A = J^^^ dre '^"''^ . The required second term in the device 



Hamiltonian Eq. ( 26 ) may now be evaluated as 



= li^^a + SV^)S^., ^ea.vep. (30) 

The approach discussed in this subsection follows the work of Pecchia et. al [37] , who 
earlier implemented the static Landauer formalism for the DFTB method. 



3.2.2. Self- energy 

From Eq. ([8|, the retarded self-energy function X)^ can be expressed in the energy 
domain as 

i:liE)=TigliE)T^, (31) 

where is the retarded Green's function of the ath contact. As in previous work, 
we assume that the contacts are semi-infinite periodic lattices, which are divided into 
principle layers (PLs) along the transport direction. The PLs are chosen so wide that 
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only interactions between nearest PLs need to be considered. Then the couphng matrix 
Tq, between contact and device region (which includes at least one PL in the device 
part) will be restricted to one PL, and only the surface block of g^^ the surface Green's 
function gl^{E)^ is needed to calculate E^. 

The surface Green's function g%{E) can be obtained from the standard renormaliza- 
tion method [38 . Denoting Hqq as the Hamiltonian of one PL, and Hqi as the coupling 
matrix between the nearest PL, g^^ is given by 

gl{E) = {EI - C^r' (32) 

where the matrices C,^ are calculated from the recurrence relations 

C+i = C^^a,{EI-r^^)-^h, (33) 

77,+! = r^,^a,{EI -r)^)-^h,^h,{EI -r),)-^a, (34) 

tti+i = ai{EI -rj^y^ai (35) 

= biiEI-rj^^bi (36) 

with the initial matrices (^q — Ha,00i'na o — Ha^oo^o^o = Ha,oiibo = -H"J^oi- 
simplest approximation, the PL matrices Ha^oo and Ha^oi are obtained from a conven- 
tional DFTB calculation in which the contact is treated as a finite cluster. Although the 
semi-infinite nature of the contacts is still accounted for by the renormalization method, 
computational artefacts have recently been reported for this approach Instead, 
we obtain the PL matrices from simulations using periodic boundary conditions with 
converged Brillouin zone sampling [40 . 

After construction of the surface Green's function gl^{E)^ it is straightforward to 
compute the self-energy at the Fermi energy 

i:i(EF)=Hl,,gliEF)H^,oi, (37) 

which provides the necessary parameters for the wide-band approximation used to de- 
termine the dissipation term Q^. 

3.2.3. Initial density matrix 

We assume that the device is in thermal equilibrium at t = with a common Fermi 
energy for the full lead-device-lead system. Only for t > a bias is applied which drives 
the central region out of equilibrium and results in a current fiow. The initial conditions 



for the solution of our central equation ( 10 ) may therefore be obtained from equilibrium 
Green's function theory. For the reduced density matrix cr(0) we have 

/ + 00 TTp 
^ —fiE)-s{G^iE)}, (38) 

where ^{G^} is the imaginary part of retarded Green's function 

G^{E) = [EI - H{0) - 5]"(^)]"^ (39) 



Notice that H{0) is determined by cr(0) through equation (26) and (30), thus the above 
equations need to be solved self-consistently. Compared with the conventional DFTB 
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scheme for closed systems, equation (38) replaces the secular equation (24) to obtain the 
density matrix, while equation (26) replaces equation ([25| to construct the Hamiltonian. 

After obtaining the initial conditions cr(0) and i^(0), we use the Runge-Kutta method 
to propagate the density matrix a{t) with a time step At of 0.02 fs. Solution of the 
Poisson equation allows the Hamiltonian H{t) to be updated according to equation 
(26). In this approach, H{t) depends only on the density at time t and not on densities 
at former times, which corresponds to the adiabatic approximation for exchange and 
correlation. Finally, the dissipation term Q^it) is calculated from equation (14) to 
compute the current and a{t-\-At). This completes the discussion on the implementation 
of the TD-DFTB method for open systems. 



4. Test calculations 



The method described above has been applied to two prototypical molecular junc- 
tions: a carbon chain with metallic transport properties and an oxygen anchored benzene 
molecule with small linear response conductance. In Sec. 4.1 we will discuss the estab- 
lishment of a steady state for asymptotically constant bias potential, while Sec. |4.2| is 
devoted to a comparison with first principles TDDFT results in order to classify the 
accuracy of our approximate scheme. 



4.1. Development of steady state 

The first system we study here is a carbon chain which bridges two aluminum wires 
as shown in Fig. [1J3. The chain consists of seven carbon atoms with a distance of 1.32 
A between nearest neighbors. The electrodes are given by Al nanowires of finite cross 
sections oriented in (001) direction. The two end carbon atoms of the chain are located at 
the hollow sites of the surfaces, and the distance between the end atom and the first layer 
of the Al surface is 1.0 A. Each PL in the electrodes comprises 18 Al atoms. Besides the 
carbon chain, two adjacent PLs are included in the device region amounting to 43 atoms 
in total. Similar systems have been studied by several authors with the DFT-NEGF 
method [4TH43J. Thus, it serves as a good test model for our purposes here. 

In an initial application of the TD-DFTB scheme, we apply a time dependent bias 
potential that is exponentially turned on with a time constant T and approaches a 
constant value Vq asymptotically: 

Vl = 0, VR{t) = V°(l-e-'/^). (40) 

Fig. [2] depicts the current traces for various time constants T. Although the initial 
transient currents differ considerably, the same final steady state current is reached in 
all cases. This is a nontrivial result, as earlier model studies indicated the possibility of 
multiple steady state solutions [44], or even the absence of a steady state [22j[45]. The 
latter investigations by Khosravi and co-workers predict persistent current oscillations if 
the device region supports a bound state which is uncoupled to the leads. 

The existence of such undamped solutions may also be deduced in the present for- 
malism. If we assume that a{t) and H{t) approach constant values at tg such that 

a{ts) ^ o-(oo) A H{ts) ^ H{oo), (41) 
11 




Figure 2: TD-DFTB transient currents through the carbon chain of Fig. ^ for an exponential turn on 
of the bias with different time constants T. The inset shows the time dependent bias approaching = 
3 V 



the last part in the dissipation term of Eq. 14 may be decomposed into 



for t > tg- Here K"^ is given by 



(42) 



2i 

TT 



de 



Ae-VO)I-Hi^)-i:^ 

and contributes to the steady state current, while the component 

with 



h.c. 



el - H{0) - S'- (e - VS) I - H{oo) - S*^ 



(43) 



- / deexp [i{{e- V°)I - H{oo) - S*^} (t - K~(e) + h.c. 

^ J —OO 



(44) 



might give rise to permanent oscillations of the current beyond t^. This occurs if at least 
one of the eigenvalues of the effective device Hamiltonian H -\- T,^ is purely real, i.e. a 
bound state exists which is completely uncoupled from the contacts. As an additional 
trivial prerequisite, al least one transmitting channel must be coupled to the leads, oth- 
erwise A and therefore also the current vanish. In all of our practical simulations so far, 
we never encountered such a situation of permanent current oscillations. 



4.2. Comparison with TDDFT 

In this section we benchmark the approximate TD-DFTB scheme against results ob- 
tained with the LODESTAR code, that implements the reduced density matrix propagation 
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Figure 3: Left figure: Time dependent currents for the carbon chain device of Fig. ^ as given by TD- 
DFTB and TDDFT using the LDA exchange-correlation functional with the two basis sets STO-3G and 
6-31G. Simulation parameters: V^^ = V, V"^ = 0.1 mV, T = 1 fs. The solid hue for the TD-DFTB 
method is the current at the left lead-device interface, while the filled triangles give the current at the 
right lead-device interface. Right figure: Transmission functions T(E) a,t V = Y obtained from the 
static DFT-NEGF approach at the respective level of theory. 



method of Sec. [2] at the fuh TDDFT level. Throughout the paper the Perdew-Burke- 
Ernzerhof exchange-correlation functional is employed for TD-DFTB and the TDDFT 
calculations are performed in the adiabatic LDA. Two different Gaussian- type basis sets 
are used in the latter: the ST0-3G set is roughly of the same size as the minimal DFTB 
basis, while the more accurate 6-31 G set is a split- valence basis of double-^ quality. 

We will first discuss the regime of linear response with small bias. Results for the 
temporal bias profile of Eq. (40) with T = 1 fs and = 0.1 mV are given in Fig. [s] 
Currents can be evaluated at either the left or right lead-device interface. In the initial 
phase these current may in principle differ in magnitude, indicating a transient charging 
and decharging of the device. In the steady state both need to be equal, which provides 
a stringent test for the numerical accuracy of the implementation. 

Comparison of the current traces shows that a steady state is reached in roughly 5 fs, 
although especially in the TDDFT/ST0-3G and TD-DFTB methods smah amplitude 
oscillations are observed which fully die out only after several tenths of fs. The fact that 
the current does not follow the bias instantaneously was related to a kinetic inductance 
in Ref. [46], which originates from the inertia or effective mass of the carriers in the leads. 
Focussing now on the long time limit and taking the 6-31 G results as reference, the steady 
state currents of TD-DFTB and TDDFT/ST0-3G show a deviation of roughly 30 %. 
The conductance of the latter is with G =1.1 Go (Go = 77.48 /iS) in good agreement 
with the single-C results of Ke et.al. [43 . 

One drawback of time-dependent transport simulations are the limited options to 
analyze the outcome. In contrast, static DFT-NEGF simulations are usually discussed 
in the framework of Landauer transport theory which features a bias and energy depen- 
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Figure 4: Left figure: Time dependent currents for the carbon chain device. Simulation parameters: Vl 
= V, = 1.0 V, T = 1 fs. Right figure: Transmission functions T(E) at V = 1 V. The vertical lines 
indicate the bias window. 



dent transmission function T(E,V) as key object [1 . The transmission provides detailed 
information on the molecular resonances, as well as on the device- lead coupling and the 
transparency of individual transport channels. A priori, it is not obvious that TDDFT 
simulations follow the Landauer picture and recent investigations show that this is indeed 
not the case [ ISHTT] . Using local or semi- local exchange-correlation functionals, however, 
steady state TDDFT and static DFT-NEGF currents agree with high precision [18 . For 
the present system we find for example a DFTB-NEGF current of 4.14 nA compared 
to a TD-DFTB value of 4.19 nA. In order to facilitate a deeper analysis. Fig. [3] also 
contains the transmission functions obtained at the different levels of theory from static 
DFT(B)-NEGF simulations in the WBA. In the linear response regime, the current is 
governed by the transmission T{Ef) at the Fermi energy, which shows the same order 
(DFT/ST0-3G > DFT/6-31G > DFTB) as found for the steady state current. The high 
transmission of the ST0-3G mainly stems from the tail of a transmission channel at 0.5 
eV above Ep, which is found at higher energies for the other two methods. To the 6-31G 
transmission at Ep contributes also a broad feature centered around -1.0 eV, which splits 
into narrower peaks at the DFTB level. This fact indicates a smaller imaginary part of 
the self-energy for this channel, which might also explain the slower decay of oscillations 
in the transient TD-DFTB current (Fig. |] left). 

Turning now to current traces at a higher bias value of F = 1 V in Fig. |4j we find 
the two TDDFT results to be close, although this is somewhat coincidental. At finite 
bias, the current in the Landauer picture is obtained by integrating the transmission 
over the bias window from V to 1 V as indicated in the figure. For the ST0-3G basis 
the resonance above Ep strongly contributes, whereas only a fraction of this channel 
resides in the bias window for DFTB and DFT in the more accurate 6-31G basis. DFTB 
resonances that are located in equilibrium at -0.25 V and -0.75 V with respect to Ep 
shift into the bias window, but have a small transmission leading to the lower value of 
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the current observed. 

The sizable variation of the transmission characteristics among the different methods 
is also seen for the next device we studied. Fig.[T]3 shows the structure of 1,4-benzenediol 
inbetween Al wires of the same kind as in the previous simulations. The DFTB optimized 
molecule was placed symmetrically into the junction so that the terminating oxygen 
atoms are positioned on the hollow sites of the Al surfaces at a distance of 1.22 A. This 
geometry is also used in the following TDDFT calculations. 

Focusing first on the left part of Fig. [5j one observes a strong initial overshoot of the 
current in all methods, before it reaches a constant value after roughly 5 fs. Again, the 
DFTB current shows longer lasting oscillations than TDDFT/6-31G, but the limiting 
values agree. Also in the transmission the DFTB results are closer to TDDFT in the 
larger basis. T(E) features a transmission gap of roughly 2.5-3 eV which is smaller than 
the corresponding gap for the isolated molecule (C6O2H6), for which we find a value of 
roughly 3.8 eV in all methods. The highest occupied molecular orbital (HOMO) is fully 
transparent at the DFTB and DFTB/6-31G level and located at -1.5 eV in Fig. [sj with 
a considerably reduced width in the former method. Inspection of the local density of 
states (LDOS, not shown here) reveals that also the DFT/ST0-3G scheme has a state 
at -1.9 eV, which carries however no current and does not appear in the transmission. 
The precise location of the lowest unoccupied molecular orbital (LUMO) is more difficult 
to access. The LDOS exhibits peaks at approximately 1.5 eV (DFTB), 0.9 eV and 2.0 
eV (DFT/ST0-3G) and 0.5 eV, 0.9 eV and 1.5 eV (DFT/6-31G), respectively. Judging 
from the LDOS alone it is difficult to discriminate between molecular states and metal 
induced gap states that appear due to the presence of the additional Al layers in the 
extended molecule. The latter are usually localized and do not significantly contribute 
to the transmission. Taking this into account, we therefore tentatively assign the LUMO 
state to the resonance at ~ 1.5 eV for DFTB and DFT/6-31G as well as ~ 2.0 eV 
for DFT/ST0-3G. These observations fit into the general trend observed for gas phase 
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molecules: The DFTB method overestimates HOMO-LUMO gaps in comparison with 
DFT calculations with converged basis sets, but to a much lesser extent than DFT 
simulations with a basis of single-^ quality [28]. The HOMO-LUMO gap has a strong 
impact on the transport characteristics of a devic^and it is known that DFT simulations 
overestimate conductances. This is due to the fact that the DFT HOMO-LUMO gap is 
significantly smaller than the true quasiparticle gap, given by the difference of ionization 
potential and electron affinity. In this regard, one would expect that the currents obtained 
from DFTB are closer to experimental values. This is certainly another case of achieving 
the right result for the wrong reason, but this point could be of interest for investigations 
targeting pragmatic solutions. 

Returning to Fig. [sj it is obvious that T{Ef) and therefore also the current is very 
sensitive to the precise peak position and broadening of the molecular states. In fact, it is 
not uncommon that transport simulations using slightly different levels of theory or even 
different implementations of the same level of theory differ by one order of magnitude 
in the predicted currents [2 . In this respect, the TD-DFTB scheme falls well into these 
general error bars. 

At the end of this section, we shortly report the required CPU time for the simulations 
presented. On a single core of an Intel Xeon X5560@2 . 80GHz multi-core processor, 
100 time steps for the carbon chain device took roughly 120 CPU minutes (TDDFT/6- 
31G), 43 minutes (TDDFT/ST0-3G) and 8 minutes (TD-DFTB), respectively. For the 
1,4-benzenediol device, the simulations took 139 minutes (TDDFT/6-31G), 48 minutes 
(TDDFT/ST0-3G) and 7 minutes (TD-DFTB), respectively. 

5. Summary and outlook 

In this article, we presented theory, implementation and validation of an approximate 
TDDFT method for open systems out of equilibrium. The temporal profile of the TD- 
DFTB current traces was found to be qualitatively similar to first principles simulations. 
A rapid switch leads to an overshoot of the current which settles into a steady state after 
only a few fs. As in earlier TDDFT simulations [18 , a steady state is always reached and 
does not depend on the history of the applied bias. Notwithstanding, we showed that 
permanent current oscillations are valid solutions of our employed equations of motion. 
The absolute value of the steady state currents was then discussed in the language and 
framework of the conventional static DFT-NEGF formalism. We found that TD-DFTB 
and TDDFT can differ significantly from each other, but this difference is not larger 
than the one between full TDDFT calculations using different basis sets. With respect 
to computational efficiency, the TD-DFTB approach is roughly one order of magnitude 
faster than a first principles implementation in the same formalism. 

A limitation of the presented approach is related to applications which require a large 
bias voltage of several V. In this case also molecular resonances far away from the frontier 
orbit als are moving into the bias window which, especially for the unoccupied states, are 
insufficiently described at a minimal basis set level. As mentioned earlier, it is entirely 



^More precisely, the equilibrium conductance is most strongly influenced by the level closest to Ep. 
Since Ep is often located in the middle of the HOMO-LUMO gap, the size of the latter is also determining 
the transport gap. 
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possible to include higher angular momentum states in the DFTB basis, although the 
two-center approximation for the Hamiltonian matrix elements becomes less accurate in 
this case. A larger basis does therefore not necessarily lead to a better result. 

The main advantage of the TD-DFTB scheme is its computational efficiency. This 
allows for important extensions which are not easily realizable in a first principles frame- 
work. To this class belongs for example the implementation of the hierarchical equa- 
tion of motion approach [24 , which goes beyond the wide band approximation for the 
leads. Conceivable is also the consideration of quasiparticle corrections to the transport 
gap. For the static case, such GW calculations have already been performed within the 
DFTB framework [W, ^T. Having applications to molecular photoswitches in mind, an- 
other important extension would be the realization of molecular dynamics simulations 
under current flow involving time-dependent external fields. Efforts in this direction are 
currently under way. 
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